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Abstract 

In this paper we reconsider the sub-Riemannian cortical model of image completion introduced in 
[16, 61]. This model combines two mechanisms, the sub-Riemannian diffusion and the concentration, 
giving rise to a diffusion driven motion by curvature. In this paper we give a formal proof of the existence 
of viscosity solutions of the sub-Riemannian motion by curvature. Furthermore we illustrate the sub- 
Riemannian finite difference scheme used to implement the model and we discuss some properties of the 
algorithm. Finally results of completion and enhancement on a number of natural images are shown and 
compared with other models. 


1 Introduction 

In this work we aim to face the problem of completion (also known as inpaiting) and contours enhancement: 
both require an analysis of elongated structures, such as level lines or contours. The word inpainting^ 
first coined by Bertalmio, Sapiro, Caselles, and Ballester in [3], refers to a technique whose purpose is to 
restore damaged portions of an image. In terms of psychology of vision inpainting corresponds to amodal 
completion^ a phenomenon described by Kanizsa in [44] , which consists in the completion of the occluded part 
of an image. This problem can be faced with many techniques such as geometric and variational instruments 
which start from the selection of existing boundaries or algorithms acting directly on the image space, based 
on the Gestalt law of good continuation and the classical findings of psychology of vision. This last class of 
algorithms goes back to the works of Kanizsa [44] , Ullman [65] , Horn [39] and find the prototype of models 
for curve completion into Mumford’s Elastica functional (see [50]): 

f \/l + ds. 

A 

This expression minimizes a functional defined on the curve itself, where k is the curvature along the curve 
7 , parametrized by arc length s. The models proposed by Mumford, Nitzberg and Shiota, see [50] and [52], 
and by Masnou and Morel, see [48], generalize this principle to the level sets of a gray-valued image. In 
fact, if is a squared defined in and / : I) M is an image, the considered functionals are second order 
functionals defined on /: 



A model expressed by a system of two equations, one responsible for boundary extraction, and one for figure 
completion was proposed by Bertalmio, Sapiro, Caselles, and Ballester in [3]. In Sarti et al [63] the role of 
the observer was considered, letting evolve by curvature in the Riemannian metric associated to the image 
a fixed surface called point of view surface. 

*Dipartimento di Matematica, Universita di Bologna, (giovanna.citti@unibo.it) 

^Center of Mathematics, CNRS-EHESS, Paris (benedetta.franceschiello@ehess.fr) 

^Department of Mathematics and Computer Science, TUE. (g.r.sanguinetti@tue.nl) 

§ Center of Mathematics, CNRS-EHESS, Paris (alessandro.sarti@ehess.fr) 


1 







Close to the completion problem there is the one of enhancement, allowing to make the structures of 
images more visible, while reducing the noise. Also for contours enhancement, models expressed as non linear 
PDE dependent on the gradient were proposed in the same years. The first models were due to Nitzberg and 
Shiota [53], Cottet and Germain [17]. Later on Weickert [68, 69] proposed a coherence-enhancing diffusion 
method, called CED. Let us also mention the models of Tschumperle [64], who included curvature into the 
diffusion process in order to improve enhancement. 

A different class of algorithms to perform both completion and enhancement is directly inspired by the 
functionality of the visual cortex, which is expressed in terms of contact structures and Lie groups with a 
sub-Riemannian metric. We refer to Section 2 for a precise definition of this metric and for references of the 
mathematical aspects of the problem. We only recall here that PDEs in this setting are totally degenerate at 
every point, so that classical results on PDE cannot be applied to study their solutions. Then we will follow 
the approach provided by Evans and Spruck in [28]. The first geometric models of the functionality of the 
visual cortex date back to the papers of Hoffmann [36], Koenderink [46], Zucker [70]. Petitot and Tondut 
in [56] proposed a model of single boundaries completion through constraint minimization in a contact 
structures, obtaining a neural counterpart of the models of Mumford, see [50]. Sarti and Citti in [61] and 
[16] proposed to model the cortical structure as Lie group SE{2) of rigid transformations of the plane with a 
sub-Riemannian metric. Each 2D retinotopic image is lifted to a surface in this higher dimensional structure 
and processed by the long range via a diffusion driven motion by curvature. This model performs amodal 
completion on the whole image and can be considered the cortical counterpart of the models of Mumford, 
Nithberg and Shiota [52]. A large class of algorithms of image processing in Lie groups were developed by M. 
Van Almsick, R. and M. Duits and B.M. Ter Haar Romeny in [24]. The main instrument he developed is an 
invertible map from the 2D retinal image to the feature space expressed in terms of Eourier instruments. He 
faced both the problem of contour completion and enhancement in presence of crossing and bifurcating lines. 
Let us also recall the results of Boscain [9] , who implemented the diffusion step of the model introduced in 
[61] using instruments of Group Eourier transform, similar to the ones proposed by Duits in [25]. 

The paper is organized as in the following. In Section 2 we reconsider the completion model of Gitti 
and Sarti (see [61]) and propose to extend its application to perform contours enhancement. Indeed we 
will see that the action of the sub-Riemannian diffusion is strongly directional, and it reduces noise only 
in the direction of contours. Vice versa it performs no diffusion in the orthogonal direction, preserving the 
structures of the images. 

Since the model is formally expressed as a diffusion driven motion by curvature, in Section 3 we consider 
the mean curvature equation in the SE{2) group and prove the existence of viscosity solutions. Existence 
results for this type of equations in the Euclidean setting are well known (see for example [28]) but in the sub- 
Riemannian setting the existence of solution was known only for the Heisenberg group (see [10], [29], [10]). 
Here we start filling the gap, extending the existence results of [28] in the present setting. In particular we 
look for sub-Riemannian viscosity solution as limit of the analogous problem in an approximating Riemannian 
setting. The main difficult we have to face is the fact that in the sub-Riemannian setting the derivatives 
are replaced by directional derivatives, which in general do not commute. Hence we will have to introduce 
a new basis of the tangent space which commutes with the natural one. 

In section 4 we discuss the numerical scheme used to implement the mean curvature flow. Einally 
in section 5 we conclude the paper showing some numerical results of completion and enhancement, and 
comparing with different methods. 

2 Sub-Riemannian operators in image processing 

In this section we first recall the cortical model of image completion proposed by Gitti and Sarti in [61] and 
[16]. By simplicity we will focus only on the image processing aspects of the problem, neglecting all the 
cortical ones. Then we show that this mechanism can be adapted to perform contour enhancement. 
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Figure 1 : A contour represented by the curve ^ 2 D{t) is lifted into the roto-translation group obtaining the 
red curve 731 ) (t). The tangent space of the roto-translation group is spanned by the vectors Xi and X 2 . 

2.1 Processing of an image in a sub-Riemannian structure 

2.1.1 Lifting of the level lines of an image 

The cortical based model of completion, proposed by Citti and Sarti in [61] and [16], lifts each level line of 
a 2 D image I{x,y) in the retinal plane to a new curve in the group SE{2) = x 5^ of rigid motions of 
the Euclidean plane, used to model the functional architecture of the primary visual cortex (VI). Precisely 
if a level line of I is represented as a parametrized 2-dimensional curve j 2 D = and the vector 

Xi = (cos(^(t)), sin(^(t))) is the unitary tangent to the curve 72 ^ at the point t, we say that 0{t) is the 
orientation of j 2 D at the point t. Then the choice of this orientation lifts the 2D curve ^ 2 D to a new curve: 

(x(t),y{t)) {x(t),y(t),0(t)) G X 5^ (1) 

as it is shown in figure 2.1.1. By construction the tangent vector to any lifted curve = {x{t)^y{t)^0(t)) 
in SE(2) can be represented as a linear combination of the vector fields: 

Xi = co^{0)dx + ^hi(0)dy 

X2 = de. 

In other words we have associated to every point of x the vector space spanned by vectors Xi and 
V 2 . 


2.1.2 The sub-Riemannian structure 

We will call horizontal plane and denote it as i^M the tangent plane generated by Xi and X 2 at every point. 
Let us also note that there are no lifted curve with a non-vanishing component in the orthogonal direction 

X 3 = — sin(6>)9a, -h co^{0)dy. 

However derivations in the direction X 3 can be recovered by means of the commutator: 


X3 = X1X2 - X2X1 = [Xi, X2] = - sin(^) 5 , + co^{e)dy. 


This condition ensures that Xi , X 2 and their commutators of any order span the Euclidean tangent space of 
SE(2) at every given point, i.e. they satisfy the Hormander condition, see [38]. Then the structure obtained 
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with this lifting process is sub-Riemannian. On the plane HM we define the metric go which makes Xi and 
X 2 orthonormal. Hence, if a vector a = aiXi + a 2 X 2 G HM^ its horizontal norm is: 

|a|o = GW + W- (2) 

The first classical properties of the distance in these spaces have been established by Nagel, Stein and 
Wainger (see [51]), and Gromov (see [33]). We refer to Hladky (see [34]) and the references therein for recent 
contributions. 

In this setting the vector fields play the same role as derivatives in the standard setting. Hence we will 
say that a function u :M? x ^ M is of class in the sub-Riemannian sense (we will denote it as G C^j^) 
if there exists Xiu and X 2 U and they are continuous. In this case we will call horizontal gradient Vq: 

Vo^ = {Xiu)Xi + {X2U)X2. 

From the definition stated before it follows that the norm of the horizontal gradient is: 

\Vou\ = v/(Xiw)2 + (X2W)2. (3) 

In other words the horizontal gradient is the projection of the standard Euclidean gradient of u on the 
horizontal plane HM. 

2.1.3 Lifting of the image to a regular surface 

Since each level line of the image / is lifted to a curve in the 3D cortical space, the whole image is lifted to 
a graph 

{x,y) {x,y,0{x,y)). 

Clearly the graph of 0 can be interpreted as the zero level set of the function u 

u{x,y,e) = e - 0{x,y), 

and it can be identified as a regular surface in the sub-Riemannian structure. The notion of regular surface 
S was first introduced by Franchi, Serapioni and Serracassano in [31]: 

s = {{a;, y, 0) : u{x, y,0) = O and Wou{x, y, 0) 7 ^ 0}. (4) 

Since the sub-Riemannian surface S is union of horizontal curves, we say that it is foliated in horizontal 
curves. The horizontal normal of S is defined as 

_ 

~ |Von| ■ 

Note that in a smooth surface there can be points where the Riemannian gradient is not 0, but its projection 
on the HM plane vanishes: 

Xqu = 0 . 

Points which have this property are called characteristics and the normal is not defined at them. However 
these points are not present in lifted surfaces. 

2.1.4 Diffusiou aud couceutratiou algorithm 

We have seen in subsection 2.1.3 how to lift an image /(x,^), to a surface S. After that let us lift the level 
lines of the image I{x^y) to the function 

v{x,y,0{x,y)) = I{x,y) 

defined on the surface. The surface S and the function v defined on S will be processed through differential 
operators defined on SE(2)^ which model the propagation of information in the cortex. More precisely two 
mechanisms operate on the lifted surface S\ 
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(a) a sub-Riemmanian diffusion along the vector fields Xi and X 2 which model the propagation of infor¬ 
mation through the cortical lateral connectivity. This operator can be expressed as 

- X^ - Xi 

where Xf expresses the second derivative in the direction Xi. The operator is formally degenerated, 
in the sense that its second fundamental form has 0 determinant at every point. It has been deeply 
studied starting from the classical works of Hormander in [37], Rothshild and Stein in [58] and Jerison 
[43] and it is known that it is hypo-elliptic. After that a large literature has been produced on these 
type of operators, and we refer to [12] for recent presentation of the state of the art. 

(b) a concentration on the surface of maxima to model the non-maximal suppression mechanism and the 
orientation tuning. 

In the Euclidean setting Merrimann, Bence and Osher, proved in [49] the convergence of a similar two step 
algorithm to the motion by curvature. In Citti and Sarti [61] and [16] the authors studied the motion when 
(a) and (b) are applied iteratively and proved that at each step the surface performs an increment in the 
normal direction with speed equal to the sub-Riemannian mean curvature. 


2.1.5 Mean curvature flow 

The notion of curvature of a surface at non characteristic points is already well understood, see ([20], 
[35], [14], [57], [12]). It can be defined either as first variation of the area functional, either as limit of the 
mean curvature of the Riemannian approximation or as horizontal divergence of the horizontal normal: 


Ko = divo(z^o) = divo 


W^\)' 


If each point of the surface evolves in the direction of the normal vector with speed equal to the mean curva¬ 
ture, we say that the surface is evolving by mean curvature. From the previously expression of the curvature 
we formally get the following equation for the flow, which we can call horizontal (or sub-Riemannian) mean 
curvature flow: 



where Sij is the Kronecker function. An existence result for this equation was not known. We will provide in 
next section an existence theorem which will allow us to handle also characteristic points, and is expressed 
in terms of viscosity solutions. 


2.1.6 Laplace-Beltrami flow 


Citti and Sarti also conjectured that as a result of the previous mechanisms the function v{x,y,0{x,y)), 
which contains the gray-levels values, evolves through the flow described by the Laplace Beltrami operator 


V (5- ■- x^X^ 

2^ |Vou|2 J^i^j 

i,j = l ^ ^ 


\ |Vo«p 

v(-,0) = Vq. 

From now on in order to simplify notations we will denote: 


V in n c X 5 ^ 


^ XfuX^u 

A°JVou) = 6ij - 1 ^ , , i,j = 1,2. 


|VoM 


(6) 

(7) 


Let us note that the described equations become degenerate and the solutions are regular only along the 
directions of the foliation. 
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2.2 Enhancement and Inpainting in Sub-Riemannian geometry 

2.2.1 Inpainting of missing parts of the image 

In the previous section we described an algorithm proposed in [61] for restoring damaged portions of an 
image, where the corrupted set uj is known a priori. 

1 An image I{x^y) are lifted to a surface S = {{x^y^0{x^y)} in the Lie group SE(2) of rotation and 
translation, and the gray level of the image I{x,y) to a function v{x,y,0{x,y)) defined on S. In the 
lifting the corrupted part of the image becomes = cj x 5'^, where no surface is defined. 

2 The surface S and v{x, y^ 0{x^ y)) are processed via the algorithm of diffusion and concentration in the 
corrupted region H, where we impose Dirichlet boundary conditions. This leads to the motion by mean 
curvature of the surface S and to a Laplace Beltrami flow for v{x,y,0{x,y)). 

3 The final result is obtained by re-projecting onto the plane of the image the value of the intensity 
v{x,y,e{x,y)). 

The algorithm has been implemented in [60] via a diffusion and concentration method, while it has been 
implemented via the curvature equation in [61]. 

2.2.2 Enhancement of boundaries 

One of the scope of this paper is to extend the previous completion algorithm to solve the problem of contours 
enhancement. The aim of this technique is to provide a regularization in the direction of the boundaries, 
making them clearer and brighter and eliminating noise. We refer to the paper of Duits [25],[26] for some 
results of image enhancement in this space. Precisely he lifts the image / on the 3D features space, using 
an invertible map defined through Fourier analysis. The lifted version of the image / is processed in the 3D 
space and then reprojected on the 2D plane to recover an enhanced version of the image /. In particular he 
also provides results of enhancement in presence of bifurcation or crossing. In this paper we face the same 
problem adapting the algorithm recalled in the previous section. 

1 First we lift the level lines of an image I{x^y) to a surface S = {{x^y^d{x^y))} and we lift the gray 
levels of I{x^y) to a function v{x,y,0{x,y)) always defined on S. 

2 Then we process the surface S via a mean curvature flow and v via a Laplace-Beltrami flow. In order 
to perform enhancement we propose here to let equations (5) and (6) evolve in the full domain M? x S^. 
Let us remark that lifting the image in the 3D group allows to solve the problem of crossing elongated 
structures. Indeed if two lines cross in the 2D space and have different orientations, they are lifted to 
the 3D space to two different planes, allowing completion and enhancement. The directional diffusion 
will give place to a regularization only in the direction of contours. 

3 Finally we project into the plane of the image the value of the gray intensity v{x,y,0{x,y)). 

3 Existence of viscosity solutions 

In this section we provide the main result of this paper which is the proof of existence of viscosity solutions 
for the mean curvature flow in SE{2). We explicitly note that we do not need to develop new results for the 
Laplace-Beltrami operator, which is linear. 

As we can immediately observe the PDF becomes degenerate in the singularities of the horizontal gradient 
of the solution Furthermore, just as in the Euclidean space, we cannot expect the smoothness of 

the solution to be preserved for all times. The notion of viscosity solution has been introduced in order to 
overcome these type of problems. The method of generalized (viscosity) solutions independently developed 
by Chen, by Giga and Goto [15], by Evans and Spruck [28], by Crandall, Ishii and Lions [18] is now applied to 
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large classes of degenerate equations [42] . Recently it has been extended also in the sub-Riemannian setting, 
see [67] [10] [29]. Finally let us recall that Dirr et ah [23] have recently studied a probabilistic approach to 
the mean curvature flow in the context of the Heisenberg group. 


3.1 The notion of viscosity solution 

3.1.1 Viscosity solutions in Jet spaces 

The cortical model previously discussed, associates to each 2D curve 72 ^ its orientation. This procedure can 
be considered as a lifting of the initial image /(x, y) to a new function u in the Jet-space x of position 
and orientation. We refer to Petitot and Tondut, who first described the analogous cortical process as a 
lifting in a jet space [56]. In this section, we will sometimes denote ^ = {x,y,0) as an element of x S^. 
It is then natural to lift the function u into another Jet-space which contains the formal analogous of its 
sub-Riemannian gradient and the formal analogous of its second derivatives X^Xj. The definition of 
viscosity solution in Jet-spaces is based on the Taylor expansion, expressed in terms of these differential 
objects. Since the analogous of the increment in the direction of the gradient p is expressed through the 
exponential map, then the increment from a point ^ in the direction Vi^i is expressed as 

2 

w(exp(^7?iX°)(e,t+ s) 

At non regular points, such as kinks, there is not either a unique vector p which identifies the horizontal 
gradient and a unique matrix Vij which identifies the horizontal Hessian. Hence we need to give a more 
general notion: a couple (p, r) where = 1,2 denotes an horizontal vector and (vij) a 2 x 2 matrix is a 
superjet for u if it satisfies the following formal analogous of the Taylor development: 

2 2^2 

u{expC^riiXi){£,),t + < '^Pi'rn + 2 + <1^ + o{\ri\^ + s^). (8) 

*=i *=i i,i=i 

Let us note that if the superject exists it can be used in place of the derivative; furthermore a function u 
is a Jet-space viscosity solution if the differential equation in which the derivatives are replaced with the 
elements of the superjet is satisfied. More precisely: 

Definition 3.1. A function u G (7(M^ x 5^ x [0, oc))n>C^(M^ x 5^ x [0, oc)) is a jet space-viscosity subsolution 
of equation (5) if for every {pn^'^ij) i'^ the super-Jet we have: 


Z)L=i ^%{ph] 


9< 




if \Ph\ y 0 

for some |p| < 1 , if \ph\ = 0 . 


(9) 


An analogous definition is provided for a viscosity supersolution. Then a viscosity solution is a function 
which is both a subsolution and a supersolution. 


3.1.2 Viscosity solutions via test functions 

The definition of viscosity solution in Jet-space of a second order equation can be identified as the approxi¬ 
mation of the solution u via a second order polynomial, whose coefficients are exactly the elements (pi^rij) 
of the Jet space. More generally we assume that the function u is locally approximate by a smooth test 
function (f. This definition imposes the behavior of the function u at points where u — <f attains a maximum. 
At such points u and (j) will have the same first derivatives, so that Vo0 results to be an exact evaluation of 
the approximation of Looking at second derivatives, it follows that for every i we have: 

XiXi{u - 4>) <0, 

SO that the curvature of (f is an upper bound for the curvature of u. Due to this observations we can give 
the following definition: 
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Definition 3.2. A function u G (7(M^ x 5^ x [0,oo)) is a viscosity subsolution of (5) in x 5^ x [0, oo) 
if for any (^,t) in x 5^ x [ 0 , 00 ) and any function (j) G C{M? x x [0,oo)) such that u — <f has a local 
maximum at (^, t) it satisfies: 


M< 


EL=i {^a4>)XiX^4>, if I Vo(/.| ^ 0 

Ei,j=i ^%ip)XiXj(t), for some p e \p\ 7 ^ 1, ij\Xo(l>\ = 0 


A function u € x 5^ x [0, 00 )) is a viscosity supersolution of (5) if: 


dt(p > 


Elj=iA%{Vo<P)XiXjcP z/|Vo(/.|7^0 

EL=i ^%{p)XiXict> for some p e \p\ 7 ^ 1 , z/|Vo<^| = 0 


( 10 ) 


( 11 ) 


Definition 3.3. A viscosity solution of (5) is a function u which is both a viscosity subsolution and a 
viscosity supersolution. 


Theorem 3.4. 


The two definitions of jet spaces viscosity solution and viscosity solution are equivalent. 


3.1.3 Vanishing viscosity solutions 

A vanishing viscosity solution is the limit of the solutions of approximating regular problems. 

Let us first explicitly note that the coefficients Aij are degenerate: when the gradient vanishes, they are 
not defined. Hence we will apply the regularization procedure proposed by Evans and Spruck in [28] to face 
singularities, which consists in replacing the coefficients with the following ones: 


AIM 



PiPj \ 
|p|2 + ry ■ 


This approximation has a clear geometric interpretation, already provided by Evans and Spruck. In equation 
(5) each level set of u evolves by mean curvature. What we obtain adding a new parameter is the evolution 
of the graph of u 

Tf = M^n+l)&RM^n+l=uiM} 

and the introduction in the space of a metric depending on r. In this approximation equation (5) reads as: 


/ ut = E AJJXou)XiXjU in ncR^ xS^ 

\ hi=i 

1 t(-,0) = Uq. 


( 12 ) 


We will now introduce a Riemannian approximation of the mean curvature flow in the graph approxima¬ 
tion we made before. We extend go on the whole space SE{2) to a metric which makes the vectors Xi, 
X2, eXo orthonormal. Let us note that g^ is the Riemannian completion of the horizontal metric. Prom now 
on, in order to simplify notations we will always denote 


X^ = Xi, X| = X2, X3^ = eX3. 


(13) 


The Riemannian gradient associated to the metric will be represented as: 


VeU = XluXl + XluX^2 + X^uXl 
and, using the fact that X^ are orthonormal, we get: 

|VeU| = V(Xiu)2 + (X2U)2+e2(X3U)2. 

In the Riemannian setting equation (12) reads as: 

/ wt = E A^X(V,u)XIX^u in fl c m2 X 51 
■j hi=i 
u(-,0) = uo 


(14) 


(15) 
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where 


A\nsi,u) 



XfuX]u \ 
|VeW|2 + rJ ■ 


In order to prove the existence of a solution we apply another regularization, always introduced by Evans 
and Spruck. It consists in adding a Laplacian, ensuring that the matrix of the coefficients have strictly 
positive smallest eigenvalue. Then the approximated coefficients will be: 


Al’p-{p) = Al';{p) + aS,, 

and the associated equation becomes: 

I ut= in 

[ 'u(-,0) = Uq. 

This condition makes the coefficients satisfy the coercivity condition and allows to apply the standard theory 
of uniformly parabolic equations. 

We are now in condition to give a third definition of vanishing viscosity solution, 

Definition 3.5. A function u is a vanishing viscosity solution of (5) if it is limit of a sequence of solutions 
y^ek,rk,(Jk of equation (16). 

We will see that any Lipschitz continuous vanishing viscosity solution is a viscosity solution of the same 
equation (see Theorem 3.12 below). 


3.2 Solution of the approximating equations 

The aim of this sub-section is to prove the existence of solutions for the approximating equation (16) and, 
even more important, to establish estimates independent of all parameters, which hold true also the limit 
equation (5). We will study solutions on the whole space, since this is the assumption needed for the 
enhancing algorithm. 

Theorem 3.6. Assume that uq G (7^(M^ x S^) and that it is constant on the exterior of a cylinder, i.e. 
there exists S > such that: 


uq is constant on x 5^ D {x^ y^ > S'}. (17) 

Then there exists a unique solution G xS^ x [0, oo)) of the initial value problem (16). Moreover, 

for all t > 0 one has: 

t)||£oo(]^2x5l) < ||'?^o||£'^(R2x5i) (18) 

||V£;'U^’^’'^(-, t)||£oo(M2x5l) < ||V£;'Uo||£oo(M2x51)- (19) 

where Ve(’) denotes the Euclidean gradient. 

Since the previous estimates do not depend on a and e, they will be stable when these parameters go to 

0 : 

Corollary 3.7. The solution M of the equation (12) satisfies conditions (18) and (19). 

This results generalizes to SE{2) the previous results of [28] and [10]. The main difficult to face in the 
extension is the fact that the vector fields does not commute, hence it is not easy to find a nice equation 
satisfied by the gradient. Following the approach proposed by Mumford in [50], we will take the derivatives 
along the direction of green a family of vector fields which are right invariant with respect to 

the group law. It is a general fact that these vector fields commute with the left invariant ones. We recall 
this result in the special case of our vector fields, for reader convenience: 
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Lemma 3.8. A right invariant basis of the tangent spaee ean he defined as follow: 


Yi = d, 

Y 2 = X 2 A- {x cos 0 Yysin 0)Xs Y {xsinO — y cos 0)Xi = dg — ydx + xdy 

YS = dy. 

We will see that the veetor fields defined in (13) eommute with {Fi}i=i,2,3- 

Proof We will calculate their Lie bracket: 

[Xi,Yi] = {cosOdx -A smOdy)dx — dx{cosOdx YsinOdy) 

= cos Odxx + sin Odyx — cos Odxx — sin Odxy = 0. 

Since the coefficients of I2 do not depend on 0 it is clear that 

[^2,^2] =0. 


Finally 


[Xs,Ys] = {sinOdx — cos 0dy)dy — dy{sm0dx — cos Ody) 

= sin Odxy — cos Odyy — sin Odxy + cos Odyy = 0 

The other combinations commute analogously. □ 

The first step of the proof of Theorem 3.6 is the existence of the function u and its bound: 

Theorem 3.9. Under the assumption of Theorem 3.6 on the initial datum, the initial value problem (16) 
has an unique solution G x 5^ x [0,oc)) sueh that 

t)||£oo(]^2x5l) < ||'^o||£oo(]^2x51) (20 ) 

Proof For a > 0, consider the problem associated to equation (16) on a cylinder 5(0, r) x [0,T], with initial 
data 

<’^’""(-,0) = uo, (21) 

and constant value on the lateral boundary of the cylinder. Note that coefficients A^Y^ satisfy the uniform 
parabolic condition: 

< Alp^{p)piPj ( 22 ) 

for each p,p G hence the theory of parabolic equations on bounded cylinders ensures that for every 
fixed value of the parameters there exists an unique smooth solution (see for example Ladyzenskaja, 

Solonnikov, Ural’tseva [47]). By maximum principle we have 

ll'^r ^’'^(•5)||£'^(R2 x 51) < ||'^o|Uo°(R2x5i)- (23) 

Letting r tend to oo, we obtain a solution defined on the whole x [0,T] such that 

I loo < ||'^o||£oo(R2x5i)- 

□ 


We can now complete the second part of the proof of Theorem 3.6 which involves the estimate of the 
gradient: 

Theorem 3.10. Under the assumption of Theorem 3.6 and 3.9, the solution of the initial value problem 
(16) satisfies 

||V£;'U^’^’'^(-, t)||£oo(]^2x5l) < ||V£;'Uo||£oo(]^2x5l)- (24) 
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Proof. From Theorem 3.9 we known that there exists an unique smooth solution of equation (16) and 

we only have to estimate its gradient. To this end, we can differentiate equation (16) along the directions 
{^i}i=i,2,3, and using Lemma 3.8, we obtain the following equation for Wi = for all i = 1,2,3, and 

for uj 4 ^ = - {yodx - xody)u^’^’^: 




e,r,a 




(25) 


The parabolic maximum principle applied to the previous equation yields: 


t)||£oo(]^2x5l) < ||^i'?^o||£«=(]^2x5l) 


(26) 


This implies that 

||(?a;l(^’'^’‘^(-,t)||£oo(R2xSl) + ||9y«^’'^’°'(-,t)||£oo(R2xSl) < C11 V^Mq 11 £«> (R2 X ) • 

We have now to establish the estimate of the derivative dg. For every fixed value of {xq, yo) we have 

\dgu’"’" {xo,yo,d)\ < max \Y 2 u’"’" - {yod^ - xody)u’"’'’\ < 

\v-Vor + \^o-x\ <'^ 


< max |F2 'Wo - {vodx - xody)uo\ < 

\y-yo\‘^ + \xQ-x\‘^<l 


< max \deUo\ 

\y-yo\‘^+\xo-x\‘^<i 


, ,2™,®^ ,^^\vo-y\\dxUo\+ niax \xo - x\\dyUQ\ < 

\y-yor+\xo-x\^<i \y-yor+\xo-x\^<i 


< ||V£;rto||£oo(]^2x5i) 


□ 


Let us conclude this section remarking that the proof of Theorem 3.6 is a direct consequence of the two 
Theorems 3.9 and 3.10. 


3.3 Existence result 

In order to extend to our setting Evans and Spruck’s argument in the proof of [28], as well as the proof of 
[10] we need to let go to 0 the three approximating parameters cr ^ 0, r ^ 0 and e ^ 0. Since the estimates 
we have established are uniform in all parameters, we immediately have the existence of a vanishing viscosity 
solution: 

Theorem 3.11. Assume that uq G (7(M^ x S^) is Lipschitz continuous and satisfies (17). Then there exists 
a vanishing viscosity solution u G of (5), which satisfies the following properties: 

||'^(•5^)||£°^(R2x5l) < C'||'^o||£^(R2x5i) (27) 

||V£;'U(-, t)||£oo(]^2x5l) < C'||V£;'Uo||£oo(]^2x5i)- (28) 

Proof. Since uq is constant at infinity, we immediately deduce from Weiestrass theorem that the Euclidean 
gradient \/eUq is bounded. Employing estimates (18),(19) and Ascoli Arzela Theorem we can extract two 
sequences {<Jk}^ {^k}^ {xk} 0 of positive numbers such that ^ ^ 0 and such that the corresponding 
solutions {u^ = convergent in the space of Lipshitz functions. Then, by definition the limit 

is a Lipshitz continuous vanishing viscosity solution. □ 

We will now prove that this vanishing viscosity solution is indeed a viscosity solution: 

Theorem 3.12. Assume thatu^ G C{M? x S^) is continuous and satisfies (17). Then the vanishing viscosity 
solution detected in Theorem 3.11 is a viscosity solution u G (7^’^ of (5). 
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Proof. In order to prove that is a viscosity solution we consider a function (p G X 5^ X [0,oo)) 

and we suppose that u — p has a strict local maximum at a point (^ 0 :^o) G x 5^ x [ 0 , oo). Since u is a 
Lipschitz continuous vanishing viscosity solution, it can be uniformly approximated by solutions (u^) of the 
approximating Riemannian problem (see also Theorem 3.12). As ^ u uniformly near (^ 05 ^ 0)5 — <p has 

a local maximum at a point with 


{f,kPk) ^ {f,oPo) as k^oo 


(29) 


Since and p are smooth, we have 


Veu'" = ^E(p , dtU^ = dtp and De{u'^ - P) < 0 at (C/c,4) 


where is the Euclidean Hessian. Thus 




(30) 


This inequality can be equivalently expressed in terms of the coefficients Al’j as follows. At the point 

dtcf> - ( 31 ) 

< dtu'^ - Aiy {V^^v!‘)XrXf («'' + </>- vl^) < 0 (32) 

If Vo^(^o,fo) 7 ^ 0, also Xo(f>{^k,tk) 7^ 0 for sufficiently large k. Then letting fc —> oo we obtain from (32): 

(33) 


dtcp < Y. (% - ^,PPt )XiXj4> at {Co, to) 


^i=i 


ivo^r 


which implies that u is a viscosity subsolution. 

If Vo0(^0 5 ^o) = 0 then we set 

k _ ^ekPi^k^ tk) 


T] = 




There exists rj such that rj^ —> rj. Note that 


k\ I _ £k\X3(l>{Ck,tk)\ 


l(r)3| = 


{£k/Tk)\X3(j){Ck,tk)\ 


"-^1 _ < V o-rv,/vi 

vwzwEikW^k ppXPXXJXPXtkX 


Since the expression vanishes as /c ^ oo we have r]s = 0. The PDE (32) now reads as: 

3 

dt<p{^k,tk) - Y - ViVj)xrx;'^(i,{Ck,tk) < 0 


so as /c ^ oo we obtain 


2 

dt(t>{Co,to) < Y “ ViVj)XiXj(j){Co,to) 

hi=i 


(34) 


concluding the proof for the case in which u — p has a local strict maximum at point (^ 0 :^o)- H u — p has 
a local maximum, but not necessarily a strict local maximum at (^ 0 ,^ 0 ), we can repeat the argument above 
replacing p{x, t) with 

^{C,t) = 4>{C,t) + 1C + {t- to)'^ 

again to obtain (33),(34). Consequently u is a weak subsolution. That u is a weak supersolution follows 
analogously. 

□ 
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From the above result we can only say that there is a subsequence of which is convergent to the 

vanishing viscosity solution u. In order to prove the uniqueness of the vanishing viscosity solution, we would 
need the sub-Riemannian analogous of estimate established by Deckelnick and Dzuik in [22]: 

Proposition 3.13. There exists a eonstant C > 0 independent of a^r and e sueh that: 

(35) 

Letting e and cr go to 0 we also get: 

\\u^ - u\\<x < Ct'^ (36) 

where is a solution of ( 12 ). 


4 Numerical scheme 


In this part we provide the numerical approximation we used to implement the sub-Riemannian motion by 
curvature which performs inpainting and enhancement. Since our scheme is directly inspired by the classical 
one of Osher and Sethian (see [54]), we will explain how to adapt the discretization to the sub-Riemannian 
setting. Let us recall here how the scheme of Osher and Sethian in [54] can be adapted to our sub-Riemannian 
setting. The mean curvature flow (12) can be explicitly written as: 

^ (X2(u)) 2 . Xn(u) + (Xi(u))2 . X22{u) - X,{u)X2{u) ■ 2 X, 2 {u) X,{u)X2{u) ■ [X,,X2]{u) 

{X,{uW + {X2{u)r+T + (Xi («))2 + (X2(u))2+r 

This equation presents two distinct terms: the first part of the flow presents second order derivatives and 
corresponds to the curvature term, the second one has only first order derivatives and correspond to the metric 
connection. The solution u{x^y,0,t) is discretized on a regular grid with points Xi = iAx^yj = jAy^Ok = 
kAO, with time discretization tg = sAt. We will denoteZ)+^I/(i, j,/c, s), D~^U{i^j,k^s), D^^U{i,j^k,s) the 
forward, backward and central difference of a discrete function U at point (i, j, /c, s) with respect to x, and use 
analogous notations for the other variables y djnO. In terms of these derivatives we will define the analogous 
differences in the direction of the vector fields Xi and X 2 . Precisely if we have discretized the direction 0 
with K points, we will denote Ok = kiilK for k = 1, • • • LC, and we will call 

2 ^+Xi jjOkD^^U{i^ j, /c, s) -1- sin OkD^'^U (i, j, k, s) 


and analogously define backward and central difference D for the vector Xi and for the vector 

X 2 . Let us adapt the scheme proposed by Osher and Sethian in [54] to our case: 


(i) the first order term Xi{u)X 2 {u) • [Xi^ X 2 ]{u) is discretized using the upwind scheme for [Xi, X 2 ] 
Taking into account the upwind scheme for the vector field X 3 , the first order term is given by: 

1 _ max(- sin OkD^^^ U,0)D-^U ^ min(- sin OkD^^^ U, 0)D^^U 

^ |_D0Xi[/|2 + \D0X,JJ\2^^ + 

max (cos OkD^^^ U, 0)D-yU + min(cos OkD^^^ U, 0 )D+^I/ 


X 3 . 

(38) 


(ii) second order derivatives are implemented as usual as ^ jjOXi jjQX 2 ^ which leads 

to second order central finite difference. We will implemented as central differences the first derivatives 
coefficients of ^ . The first derivative with respect to Xi, coefficient of the 

second mixed derivative, will be upwinded as before. Generalizing an idea of [13], the denominator 
will be a mean of central derivatives: 

\DintU\\i,j, k,s)+T = ^ Y. \D°^^U\Hi,j, h, s) + ^ k, s) + r, 

ki£{k — l,k,k-\-l} ii£l 


where / is the family of indices I = {{i — 1, j), (i, j), {i + 1, j), (i, j — 1), {i^j + !)}• 
The second order discretized operator will be denoted W‘^{U){i, k, s). 
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The difference equation associated to the continuous equation (37) will be expressed as: 

k,s + l) = U{i,j, k, s) + At{W^U){i,j, k, s) + At{W^U){i,j, k, s) 

with initial condition = Uq. We recall that convergence of difference schemes for the mean curvature 

ow inspired by the scheme of Osher and Sethian has been object of a large number of papers in the Euclidean 
setting. The stability of one of them was proved in [13]. Another monotone scheme was proposed by Crandall 
and Lions (see [19]) and its convergence was proved by Deckelnick in [21] and Deckelnick & Dzuik in [22]. 

The ideas at the basis of the stability proof of [13] can be extended to the present version of the Osher 
and Sethian scheme, leading to the following result: 

,2 

Theorem 4.1. The difference scheme presented above is stable in the sense that if At < then 

||C/||oo < WUoWoo 


Proof. If [/ is a solution of the discrete equation, also V = f/ — ||[/o||oo is a solution of the same equation: 

V{i,j, k,s + l) = V{i,j, k, s) + At{W^V){i,j, k, s) + At{WW){i,j, k, s). 

Hence E(0) < 0, and we have to prove E < 0, for all time. In order to study the term we have to 

discuss the sign of ai = — sinOkD^^^V and 02 = cosOkD^^^VD^^^V: we will assume that they are 
both positive since the proof is similar in all the other cases: In this case 


{W^V)(i,j,k,s) 


ai{V{i,j,k,s) -V{i- I, j,k, s)) + a 2 {V{i, j,k, s) - l,k,s) 

\D0XiV\2 + \D0^^V\^+t 


{cos{dk)-sm{dk))DO^^V{i,j,k,s)D°^^V{i,j,k,s) V{i,j,k,s) 

{\D0^iV\^ + \DOX 2 V \2 + T)h^ 2/i2 

Analogously, having upwinded the coefficient of (IT^E)(i, j,/c, s), e get a similar behavior. The mixed 
derivatives term can be estimated as: 

2 cos( 9 /,T)O^TT>o^ 2 )y 2 ^oJV 2 y _ 10V(iJ,k,s) 

- lD.„,VP + r - ^ -P-■ 


In conclusion 

V{i,j,k,s + 1) < V{i,j,k,s){l - < 0. 

The assertion then follows by induction. □ 

Now we recall that the equation is uniformly parabolic in sub-elliptic sense. Arguing as in [22] the 
estimates of 4th order derivatives can be reduced to the estimates of graphs over the considered group. 
Hence estimates can be obtained by a recent result of Capogna, Citti and Manfredini (see [11]). Since for r 
fixed the equation is uniformly parabolic in sub-elliptic sense, these estimates, allow to prove that 

Theorem 4.2. If is the solution of (12) with initial condition uq and U is the solution of the discrete 
scheme considered here, and a is fixed, there exist a constant C = C{r, h,a) such that if At < C{t, h), then 

\u^{iAx,jAy,kAe,sAf) - U{i,j,k,s)\ < r“, 

As a consequence, applying the uniqueness Theorem 3.13, we deduce the following convergence result for 
the solution of the mean curvature equation (5) with initial condition uq 

\u{iAx, jAy, kAO, sAt) — U{i, j, k,s)\ < 


as At < C{t, h). 
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Figure 2: An example of completion performed by the algorithm. In this artificial image the image gradient 
is lifted in the x space and the black hole is completed by mean curvature flow. Since the level lines 
of the image are approximately circular, the algorithm performs very well. 



Figure 3: Completion result on a real image through sub-Riemannian mean curvature flow in X S\ as 
described in the paper. 

5 Results 

In this section we compare the results of the model with previous works, and discuss our results. The 
completion method, first proposed in [16], has now a solid theoretical background, and is compared with 
more recent results obtained by different authors. Results of our new model enhancement are proposed, 
comparing with previous results of Duits [26]. Finally we show an example of inpainting and enhancement 
of an image. 

5.1 Inpainting results 

Since the lifting procedure is based on the selection of the orientation of level lines at every point, the 
algorithm performs particularity well for completing gray level images which have non vanishing gradient at 
every point. Hence we will start with a simple artificial image of this type. The algorithm performs very 
well for completion of curved level lines. 

We will now test the algorithm on natural images. In the first image a black hole is present, and the 
algorithm correctly reconstructs the missed part of the image: 

The model ([61]) studied here performs completion via the curvature flow. Very recently Boscain et al. 
in [9], tried to replace this non linear equation by simple diffusion. In figure 4 left we consider an image 
courtesy of U. Boscain [9], partially occluded by a grid and show the results of completion performed in [9] 
(second image from left), by the heat equation on the 2D space (third image from left) and by the Citti 
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Figure 4: Left: an occluded image (courtesy of U. Boscain ([9])), second image from left: the image processed 
in ([9]); third image from left: the same image processed through the heat equation right: image the inpainted 
using the Citti Sarti algorithm. 



Figure 5: A detail of previous image: Left: the original image (courtesy of U. Boscain ([9])); Second image 
from left: the image processed in ([9]); Third image from left: the image processed through the heat equation; 
Right: image inpainted using the proposed algorithm. 


and Sarti model (right). A detail is shown in figure 5. Since the considered image is a painting, extremely 
smooth, with low contrast, the 2D heat equation is already able to perform a simple version of completion. 
In this case the implementation of sub-Riemannian diffusion [9] provides a worse result , while the curvature 
model reconstructs correctly the missed contours and level lines. 

In figure 6 (and in the detail taken from it in figure 7) we consider an other example taken from the 
same paper. In this image the grid of points which are missed is larger, and the previous effect is even more 
evident. 

In a more recent paper Boscain and al. introduced a linear diffusion with coefficients depending on the 
gradient of the initial image (see [7]), which they call heuristic. In figure 5.1 we compare the results obtained 
with this model, the heat equation on the image plane and the strongly geometric model of Sarti and Citti. 

Then we test our implementation on piecewise constant images. Since the gradient is 0 in large part 
of the image, the lifted gradient is not defined in largest part of the image. On the other side, since the 
lifting mimics the behavior of the simple cells of the VI cortical layer, the Citti and Sarti algorithm is always 
applied on a smoothed version of the image. We have applied it on a classical toy problems proposed for 
example in Bertalmio, Sapiro, Caselles and Ballester in [3]. Results are shown in figure 9. 

In figure 10 we test our method on an image taken from the survey [4]. The present reconstruction is 
correct in the part of the image characterized by strong boundaries, but the results of [4] obtained with the 
model of Morel and Masnou (see [48]) seems to be better. The main point is the boundary detection, which 
is very accurate in the model of Morel and Masnou, while here the boundaries are detected with a gradient, 
after smoothing the image. 
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Figure 6: Top left: the original image (courtesy of U. Boscain ([9])); Top right: the image processed in 
([9]); Bottom left: the image processed through the heat equation; Bottom right: image inpainted using the 
proposed algorithm. 



Figure 7: Detail of the previous image.The result of our algorithm (4th right) preserves the circular level 
lines of the iris, while simple subRiemannian diffusion (3rd right) destroy it. Simple euclidean diffusion (3rd 
right) performs at intermediate level. 



Figure 8: On the left the occluded image. From left to right: results from [7], with 2D heat equation and 
our model. 
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Figure 9: Inpainting a constant coefficient image with the Sarti Citti algorithm. 



Figure 10: On the left the occluded image. From left to right: results from [4] with the model of [48], and 
with our model. 
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Figure 11: From left to right: the original image, courtesy of R. Duits ([26]), the enhanced image using 
CED-OS, see [26] and the enhanced image obtained using the proposed method. 
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Figure 12: Left: the original image (courtesy of U. Boscain ([9])); center: image inpainted using the proposed 
algorithm; right: image inpainted and enhanced with this algorithm. 

5.2 Enhancement results 

We will show in this section results of the application of the enhancement method we have introduced in 
Section 2.2.2. Let’s recall that enhancement consists in an image filtering that underlines directional coherent 
structures. With respect to the completion problem there is no part of the image to be disoccluded and all 
the parts of the initial data are evolved. 

In Figure 11 it is shown a medical image of blood vessels to be filtered to reconstruct the fragmented 
vessels (courtesy of R. Duits ([26]). The second image from left shows the enhancement computed by using 
CED-OS, see [26], while the third image shows the result obtained using the proposed method. 

Here we finally show an example of combination of the techniques of completion and enhancement. We 
see in this case that enhancement homogenizes the original non occluded part with the reconstructed one. 
Here we propose a detail of the previous image in order to underline the effects of the discussed techniques. 

6 Conclusions 

In this paper we have proved existence of viscosity solutions of the mean curvature flow PDE in x 
with a sub-Riemannian metric. The flow has been approximated with the Osher and Sethian technique 
and a sketch of the proof of convergence of the numerical scheme is provided. Results of completion and 
enhancing are obtained both on artificial and natural images. We also provide comparisons with others 
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Figure 13: Top left: a detail of the original image (courtesy of U. Boscain ([9])); Bottom left: a detail of 
the image inpainted using the proposed algorithm; Bottom right: same detail of the image inpainted and 
enhanced with this algorithm. 


existing algorithms. In particular we have illustrated how the method can be used to perform enhancement 
and how it leads to results comparable with the classical ones of Bertalmio, Sapiro, Caselles and Ballester in 

[3] , of Masnou and Morel in [48] . In the case of image completion we compared the technique with the recent 
results shown of Boscain, Chertovskih, Gauthier, Remizov in [9]. Furthermore the method can be applied 
not only to inpainting problems but also in presence of crossing edges, hence a comparison with the results 
of edges which cannot be done using the method proposed by Bertalmio, Sapiro, Caselles and Ballester in 
[3] is now possible. 
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